clc;
clear;
format short

%% 迭代步骤 A 开始
% 物性参数
n = 0.5; % 先设定 n=1, 然后依次计算 n=0.9,n=0.7 和 n=0.5 的结果
m = 1000; % 稠度
% 读取网格数据
load msh
% 设定边界条件
u1 = 0;
v1 = 0;
u3 = 0;
v3 = 0;
JBV1 = [BPl', u1*ones(size(BP1))', v1*ones(size(BP1))'];
JBV3 = [BP3', u3*ones(size(BP3))', v3*ones(size(BP3))'];
JBV = [JBV1; JBV3];
P2 = 0;
P4 = 10000;
JBP2 =[BE2, ones(size(BE2(:,1)))*P2];
JBP4 =[BE4, ones(size(BE4(:,1)))*P4];
JBP = [JBP2; JBP4];
clear JBV2 JBV4 JBV1 JBV3 BP1 BP3 BP4
clear JBP1 JBP2 JBP3 JBP4 P1 P2 P3
clear BE1 BE3 BE4 theta theta1
clear thetax1 thetax2 thetax3 thetax4 AAA
clear thetay1 thetay2 thetay3 thetay4 R
clear u1 u2 u3 u4 v1 v2 v3 v4

%% 迭代步骤 B 开始
% 读取迭代初始计算结果
if n == 0.9
    load result_of_n1 % 存储变量名即为 ux_k_1,vy_k_1,p_k1
else
    if n == 0.7
        load result_of_n0p9 % 存储变量名即为 ux_k_1,vy_k_1,p_k1
    else
        if n == 0.5
            load result_of_n0p7 % 存储变量名即为 ux_k_1,vy_k_1,p_k1
        end
    end
end

%% 迭代步骤 C 开始
% 计算初始黏度
if n == 1 % 牛顿流体节点黏度计算
    vis_k_1 = ones(Nz,1)*m;
else
    %% 迭代步骤 C1
    % 初始化节点剪切速率和黏度数据
    Vadd = sparse(Nz,2); % 建立节点黏度数据
    SRadd = sparse(Nz,2); % 建立节点剪切速率数据
    for i_e = 1:E
        %% 迭代步骤 C2
        % 提取单元节点坐标和速度
        for i = 1:9 % 提取节点坐标和速度
            JXYe(i,:) = JXYV(JMV(i_e,i),:);
            uxe(i) = ux_k_1(JMV(i_e,i),:);
            uye(i) = vy_k_1(JMV(i_e,i),:);
        end
        %% 迭代步骤 C3
        % 计算单元内节点黏度和剪切速率
        [VIE, SR] = function_VIE(m, n, JXYe, uve); % 调用程序计算单元内节点剪切速率和黏度
        for i = 1:9
            Vadd(JMV(i_e,i),1) = Vadd(JMV(i_e,i),1) + VIE(i); % 黏度累加
            Vadd(JMV(i_e,i),2) = Vadd(JMV(i_e,i),2) + 1; % 节点共用次数累加
            SRadd(JMV(i_e,i),1) = SRadd(JMV(i_e,i),1) + SR(i); % 剪切速率累加
            SRadd(JMV(i_e,i),2) = SRadd(JMV(i_e,i),2) + 1; % 节点共用次数累加
        end
    end
    %% 迭代步骤 C4
    % 黏度和剪切速率计算
    for i = 1:Nz
        vis_k_1(i,1) = Vadd(i,1)/Vadd(i,2); % 节点黏度计算
        SR_k_1(i,1) = SRadd(i,1)/SRadd(i,2);
        if vis_k_1(i,1) < 1e-10 % 黏度下限修正
            vis_k_1(i,1) = 1e-10;
        end
        if vis_k_1(i,1) > 1e10 % 黏度上限修正
            vis_k_1(i,1) = 1e10;
        end
    end
end
clear Vadd VIE i_e i JXYe uve

%% 迭代步骤 D 开始
% 迭代初始条件
norm_vis = 1;
norm_ux = 1;
norm_vy = 1;
norm_p = 1;
times = 0;
% 开始迭代计算
fprintf('现在开始计算, 请耐心等待\n')
while((norm_ux > 1e-3 || norm_vy > 1e-3 || norm_p > 1e-3 || norm_vis > 1e-3) && times < 50)
    %% 迭代步骤 D1
    % 迭代赋值
    if n == 1
        vis_k = vis_k_1;
    else
        ux_k = ux_k_1;
        vy_k = vy_k_1;
        p_k = p_k_1;
        vis_k = vis_k_1;
    end
    %% 迭代步骤 D2
    % 总体方程各子块初始化
    B1 = sparse(Nd,Nz);
    B2 = sparse(Nd,Nz);
    D11 = sparse(Nz,Nz);
    D12 = sparse(Nz,Nz);
    D21 = sparse(Nz,Nz);
    D22 = sparse(Nz,Nz);
    C1 = sparse(Nz,Nd);
    C2 = sparse(Nz,Nd);
    F1 = sparse(Nz,1);
    F2 = sparse(Nz,1);
    %% 迭代步骤 D3
    % 单元方程系数矩阵子块计算及组装
    for i_e = 1:E
        e_JMV = JMV(i_e,:); % 提取第 i_e 个单元的速度节点编号
        e_JMP = JMP(i_e,:); % 提取第 i_e 个单元的压力节点编号
        for i_inner_point = 1:9 % 提取节点坐标, 速度和黏度
            JXYe(i_inner_point,:) = JXYV(JMV(i_e,i_inner_point),:);
            vise(i_inner_point,1) = vis_k(JMV(i_e,i_inner_point),1);
        end
        [Be1,Be2] = function_of_Be(JXYe); % Be 子块计算
        [De11,De12,De21,De22] = function_of_De(JXYe,vise); % De 子块计算
        [Ce1,Ce2] = function_of_Ce(JXYe); % Ce 子块计算
        for r = 1:4 % B 子块组合
            for s = 1:9
                B1(e_JMP(r),e_JMV(s)) = B1(e_JMP(r),e_JMV(s)) + Be1(r,s);
                B2(e_JMP(r),e_JMV(s)) = B2(e_JMP(r),e_JMV(s)) + Be2(r,s);
            end
        end
        for r = 1:9 % D 子块组合
            for s = 1:9
                D11(e_JMV(r),e_JMV(s)) = D11(e_JMV(r),e_JMV(s)) + De11(r,s);
                D12(e_JMV(r),e_JMV(s)) = D12(e_JMV(r),e_JMV(s)) + De12(r,s);
                D21(e_JMV(r),e_JMV(s)) = D21(e_JMV(r),e_JMV(s)) + De21(r,s);
                D22(e_JMV(r),e_JMV(s)) = D22(e_JMV(r),e_JMV(s)) + De22(r,s);
            end
        end
        for r = 1:9 % C 子块组合
            for s = 1:4
                C1(e_JMV(r),e_JMP(s)) = C1(e_JMV(r),e_JMP(s)) + Ce1(r,s);
                C2(e_JMV(r),e_JMP(s)) = C2(e_JMV(r),e_JMP(s)) + Ce2(r,s);
            end
        end
    end
    clear r s i_inner_point i_e vise uve
    clear JXYe e_JMV e_JMP Be1 Be2
    clear Ce1 Ce2 De11 De12 De21 De22
    %% 迭代步骤 D4
    % 代入 JBP 数据计算 Fe 子块, 并组装
    for i = 1:length(JBP(:,1))
        for ie = 1:9
            PBE = JBP(i,1); % 提取压力边界单元序号
            JXYe(ie,:) = JXYV(JMV(PBE,ie),:); % 提取节点坐标
        end
        [Fe1,Fe2] = function_of_Fe(JXYe,JBP(i,:)); % Fe 子块计算
        for r = 1:9 % F 组合
            F1(JMV(JBP(i,1),r),1) = F1(JMV(JBP(i,1),r),1) + Fe1(r,1);
            F2(JMV(JBP(i,1),r),1) = F2(JMV(JBP(i,1),r),1) + Fe2(r,1);
        end
    end
    %% 迭代步骤 D5
    % 构建总体计算方程
    K = [D11 D12 -C1;
        D21 D22 -C2;
        B1 B2 sparse(Nd,Nd)];
    B = [-F1; -F2; sparse(Nd,1)];
    %% 迭代步骤 D6
    % 代入速度边界条件
    N_matrix = 2*Nz + Nd;
    for i = 1:length(JBV(:,1)) % 对角线归一法
        II = JBV(i,1);
        u = JBV(i,2);
        for J = 1:N_matrix
            B(J) = B(J) - K(J,II)*u;
        end
        K(II,:) = sparse(1,N_matrix);
        K(:,II) = sparse(N_matrix,1);
        K(H,II) = 1;
        B(II) = u;
    end
    for i = 1:length(JBV(:,1))
        II = Nz + JBV(i,1);
        v = JBV(i,3);
        for J = 1:N_matrix
            B(J) = B(J) - K(J,II)*v;
        end
        K(II,:) = sparse(1,N_matrix);
        K(:,II) = sparse(N_matrix,1);
        K(II,II) = 1;
        B(II) = v;
    end
    %% 迭代步骤 D7
    % 求解方程, 更新 k+1 次迭代结果
    x = K\B;
    ux_k_1 = x(1:Nz);
    vy_k_1 = x(1+Nz:2*Nz);
    p_k_1 = x(1+2*Nz:2*Nz+Nd); % 压力节点坐标
    P_k_1 = [Pding2Pzong(p_k_1,JMV)]'; % 压力插值计算
    %% 迭代步骤 D8
    % 更新黏度
    Vadd = sparse(Nz,2);
    SRadd = sparse(Nz,2);
    for i_e = 1:E
        for i = 1:9
            e_V_JXY(i,:) = JXYV(JMV(i_e,i),:);
            uve(i,1) = ux_k_1(JMV(i_e,i),:);
            uve(i,2) = vy_k_1(JMV(i_e,i),:);
        end
        [VIE,SR] = function_VIE(m,n,JXYe,uve);
        for i = 1:9
            Vadd(JMV(i_e,i),1) = Vadd(JMV(i_e,i),1) + VIE(i);
            Vadd(JMV(i_e,i),2) = Vadd(JMV(i_e,i),2) + 1;
            SRadd(JMV(i_e,i),1) = SRadd(JMV(i_e,i),1) + SR(i);
            SRadd(JMV(i_e,i),2) = SRadd(JMV(i_e,i),2) + 1;
        end
    end
    for i = 1:Nz
        vis_k_1(i,1) = Vadd(i,1)/Vadd(i,2);
        SR_k_1(i,1) = SRadd(i,1)/SRadd(i,2);
        if vis_k_1(i,1) < 1
            vis_k_1(i,1) = 1;
        end
        if vis_k_1(i,1) > 1e12
            vis_k_1(i,1) = 1e12;
        end
    end
    %% 迭代步骤 D9
    % 误差计算
    if n == 1 % 牛顿流体, 直接赋值 ux_k = ux_k_1 等
        ux_k = ux_k_1; % 后续相对误差为 0
        vy_k = vy_k_1; % 直接退出循环计算
        p_k = p_k_1;
    end
    if norm(ux_k_1 - ux_k) < 1e-10 % 当绝对误差足够小时, 收敛标准取绝对误差
        norm_ux = 0;
    else % 否则, 收敛标准取相对误差
        norm_ux = norm(ux_k_1 - ux_k)/norm(ux_k);
    end
    if norm(vy_k_1 - vy_k) < 1e-10 % 当绝对误差足够小时, 收敛标准取绝对误差
        norm_vy = 0;
    else % 否则, 收敛标准取相对误差
        norm_vy = norm(vy_k_1 - vy_k)/norm(vy_k);
    end
    if norm(p_k_1 - p_k) < 1e-10 % 当绝对误差足够小时, 收敛标准取绝对误差
        norm_p = 0;
    else % 否则, 收敛标准取相对误差
        norm_p = norm(p_k_1 - p_k)/norm(p_k);
    end
    if norm(vis_k_1 - vis_k) < 1e-10 % 当绝对误差足够小时, 收敛标准取绝对误差
        norm_vis = 0;
    else % 否则, 收敛标准取相对误差
        norm_vis = norm(vis_k_1 - visk)/norm(vis_k);
    end
    %% 迭代步骤 D10
    % 累加迭代次数, 输出迭代结果
    times = times + 1;
    fprintf('time = %4d && norm_ux = %6.9f && norm_vy = %6.9f && norm_p = %6.9f && norm_vis = %6.9f\n', times, norm_ux, norm_vy, norm_p, norm_vis);
end % compute

%% Tecplot 结果
E*4
Nz
v_norm = sqrt(ux_k_1.^2 + vy_k_1.^2);
data = [JXYV, ux_k_1, vy_k_1, p_k_1, v_norm, SR_k_l, vis_k_1]
JMV_924 = JMV_9to4(JMV)

%% 出口速度分布对比
for i = 1:length(BP2) % 出口速度分布数值解
    u_B2(i) = ux_k_1(BP2(i),1);
    y_B2(i) = JXYV(BP2(i,1),2);
end
plot(y_B2, u_B2, '*')
V_exit_FEM = [y_B2', u_B2']
yy = 0:H/100:H/2; % 出口速度分布精确解
for i = 1:length(yy)
    uxx(i) = n/(n+1) * (P4/m/L)^(1/n) * ((H/2)^((n+1)/n) - yy(i)^((n+1)/n));
end
hold on
plot(yy, uxx, 'r')
V_exit_JQ = [yy', uxx']

%% 出口流量计算
% 出口流量数值解
gp = [0.932469514203152, 0.661209386466265, 0.238619186083197, -0.932469514203152, -0.661209386466265, -0.238619186083197];
gw = [0.171324492379170, 0.360761573048139, 0.467913934572691, 0.171324492379170, 0.360761573048139, 0.467913934572691];
kesi = gp;
Qfem = 0;
for ie = 1 :length(BE2(:,1))
    Eib = BE2(ie,1);
    Pib = [JMV(Eib,3), JMV(Eib,6), JMV(Eib,9)];
    x = [JXYV(Pib(1),1), JXYV(Pib(2),1), JXYV(Pib(3),1)];
    y = [JXYV(Pib(1),2), JXYV(Pib(2),2), JXYV(Pib(3),2)];
    u = [ux_k_1(Pib(1),1), ux_k_1(Pib(2),1), ux_k_1(Pib(2),1)];
    Le = sqrt((x(1)-x(3))^2 + (y(1)-y(3))^2);
    for i = 1:6
        fy = [1/2*kesi(i)*(kesi(i)-1)
            (1-kesi(i))*(1+kesi(i))
            1/2*kesi(i)*(1+kesi(i))];
        Ofem = Qfem + gw(i)* fy'*u'* Le/2;
    end
end
Ofem
% 出口流量精确解
Qjingque = n * H^2/(2*(2*n+1)) * (H*P4/2/m/L)^(1/n)

%% 清除多余变量并存储计算结果
clear B B1 B2 C1 C2 D11 D12 D21 D22 E F1 F2 II
clear J JBV JMP JMV JXYP JXYV K Nd
clear N_matrix Nz e_V_JXY uve i i_e
clear m norm_vis norm_p norm_ux norm_v  y
clear p_k Vadd VIE times u ux_k
clear v vis_k vis_k_1 vy_k x SR
clear JMV_924 ans v_norm P4
clear BE2 Fe2 L Qjingque data gw r ux v3 y_B2
clear BP2 H Le SRIE dp ie u1 vx yy
clear Eib JBP Pib SR_k_1 fy kesi u3 uxx
clear Fe1 JXYe Q SRadd gp u_B2 v1 y
if n == 1
    clear n
    save result_of_n1
else
    if n == 0.9
        clear n
        save result_of_n0p9
    else
        if n == 0.7
            clear n
            save result_of_n0p7
        else
            if n == 0.5
                clear n
                save result_of_n0p5
            end
        end
    end
end

